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Abstract 

We present an investigation of the 2D attractive Hubbard model, con- 
sidered as an effective model relevant to superconductivity in strongly in- 
teracting electron systems. We use both hybrid Monte Carlo simulations 
and existing hopping parameter expansions to explore the low temperature 
domain. The increase of the static S-wave pair correlation with decreasing 
temperature is analyzed in terms of an expected Kosterlitz-Thouless super- 
conducting transition. The evidence for this transition is weak: If it exists, 
its temperature is very low and depends weakly on the band filling near 
half filling. The number of unpaired electrons remains nearly constant with 
temperature at fixed attractive potential strength. In contrast, the static 
magnetic susceptibility decreases fast with temperature, and cannot be re- 
lated only to pair formation. We introduce a method by which the Pade 
approximants of the existing series for the susceptibility give sensible results 
down to rather low temperature region, as shown by comparison with our 
numerical data. 

PACS 71.10.Fd 74.25.Dw 75.40.Mg 

October 1996 Submitted for publication to J. Phys. I (France) 

BI-TP 96/05 
T96/ 114 



^lacaze@spht. saclay.cea.fr 
^morel@spht. saclay.cea.fr 
^bengt@physik.uni-bielefeld.cle 
'*joergs@physik.uni-bielefeld.de 



1. Introduction 



The study of models of strongly interacting electrons is very important for 
the understanding of high temperature superconductivity. The fundamental 
mechanism of this phenomenon has not yet been clearly identified, and it may 
be interesting to study effective models, which have strong local attraction 
between the electrons. The two dimensional attractive Hubbard model is 
an example of such a model whose phase diagram as well as physical proper- 
ties in the normal non superconducting phase can be compared to properties 
of actual superconducting materials of the high superconductivity class. 

The two dimensional attractive Hubbard model is a conceptually simple 
model, which at low temperatures is expected to have a Kosterlitz-Thouless 
transition into a superconducting phase, away from half filling. At half 
filling the model has further symmetries, which prevents such a transition. 
In the absence of magnetic field, the properties of the model depend on two 
parameters (apart from the temperature) namely the coupling constant of 
the attractive local interaction U and the chemical potential /i, measured 
in units of the coefficient k of the hopping term. The model is solvable in 
the free case U = and the atomic limit k = 0. In the interesting case 
where both U and k are different from zero the model can only be studied 
through developments around the two solvable cases, high temperature series 
and numerical simulation techniques. In contrast to the repulsive Hubbard 
model, with the opposite sign of U, the fermionic determinant is positive also 
for chemical potential different from zero. Thus the numerical simulations 
do not suffer from the sign problem in this case. In fact, the two models are 
related by a change of sign in U and the exchange of chemical potential and 
magnetic field. 

In this paper, using a hybrid Monte Carlo algorithm, we perform a de- 
tailed numerical study of the model in a wide range of values of the parame- 
ters. The purpose is first to examine in detail the evidence for a Kosterlitz- 
Thouless transition and, if it takes place, to estimate the critical temperature. 
This can be done by studying the s-wave pair-field correlation function, which 
should diverge at the transition on an infinite lattice. Another interesting 
quantity is the magnetic susceptibility, sensitive to the presence of single 
electrons as contrast to those bound in pairs. One may expect the suscepti- 
bility to disappear in the strong coupling limit, where more and more pairs 
are formed. However, at fixed U the temperature behavior of the probability 
that a site has a single electron and the behavior of the susceptibility are 
not the same. Therefore the behavior in temperature of the susceptibility of 
the remnant unpaired electrons is interesting to investigate. For this latter 
quantity we make use of the series expansion in k given in Ref . , and with 
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the help of the Pade approximants method, extrapolate the series to low tem- 
peratures. We compare the results with our numerical data and show that, 
at least at small and intermediate couplings, rather low temperature may be 
reached analytically. If true also in the repulsive case of the Hubbard model, 
where numerical simulations suffer from the sign problem, series expansion 
might be interesting to use. Unfortunately, no such series are available for 
the pair field correlator. 

There have been some earlier numerical investigations of this model ^ ^, 
leading in particular to estimates of the transition temperature. We com- 
pare our estimates to theirs, when there is overlap in the parameters. Our 
estimate of the transition temperature is considerably lower than in these 
investigations. In particular, for the range of temperatures we simulated, we 
do not see a really significant difference between the half filled case, where 
no transition is expected, and away from half-filling. 

In Section 2 we define the model and the observables, which we measure. 
Section 3 is devoted to the formulation of the path integral formalism needed 
for the numerical simulation. In Section 4 we discuss the algorithm employed, 
which we have chosen as the Hybrid Monte Carlo algorithm. In Section 5 we 
discuss our data on the pair field correlation, which is a direct indicator of a 
transition into a superconducting phase. In Section 6 we present results for 
the susceptibility, comparing our numerical results with analytic results ex- 
trapolated from the series expansion. Section 7, finally, contains a summary 
and our conclusions. 



2. The Model 

The model is defined by the Hubbard Hamiltonian 

(x,y) ^ 

-/^ + , (1) 

X 

where x and y are sites on a two dimensional A'^s x Ng square lattice with V = 
sites, Qx and are coordinate space annihilation operators for spin-up 
and spin-down electrons respectively, and k is the nearest neighbor hopping 
parameter. The coupling U > denotes the strength of the attractive local 
interaction and /i is the chemical potential, defined so that /i = at half 
filling, i.e. where the total particle number (A^) = V. Here 

N = J2iala., + b%) . (2) 
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We also define the particle number density operator 

n = ^N . (3) 

We do not couple the system to an external magnetic field. Thus the 
model is invariant under an SU{2)spin- For = it is also invariant under 
another SU{2) group 0, generated by the operators 

X 

J', = \{N-V) . (4) 
The thermodynamics of the model is given by the partition function 

Z = Tt (e-'^^) . (5) 

In the simulation we will measure average values of operators in this 
ensemble, 

(O) = Tr(Oe-^^)/Z . (6) 

The basic equal time correlation functions, which give information about 
the properties of the model are the S-wave on site pairing correlation function 

P{x-y) = {ialbl + bA)ialbl + byay)) , (7) 
the magnetization density correlation function, or magnetic susceptibility 

X{x -y) = (3(^{ala^-blb^){alay-blby)) , (8) 
and the charge density correlation function 

C{x -y) = ((4a, + b%){alay + fetfo,)) . (9) 
The probability that a site is singly occupied is 

S, = ^ . (10) 

Furthermore 

P(0) = l-^i (11) 
is the probability of zero or double occupancy. 



The correlation function in Fourier space of the quantity E denoted by 
E is defined as 

E{q)=Y.e"''Eiz) . (12) 

z 

An indicator of the diverging correlation length at the phase transition, 
which has been used e.g. in 0, |^ is 

Po = P{q={0,0)) . (13) 

In the same way, the charge density wave indicator is C{tt, tt) and the uniform 
spin susceptibility is given by x(0, 0). 

From the SU{2y symmetry, it follows that at = the Hamiltonian is 
invariant in particular under the transformation 

{a^,b^) ^ {a^,P^) (14) 

with 




Using that /i = implies half filling, i.e. (n) = 1, this leads to 

P{z) = {-r[C{z)-l] . (17) 

The model is soluble in the two limiting cases U = (the free case) and 
k = (the atomic limit). For reference, these solutions and the corresponding 
values of the correlation functions are given in Appendix A. 



3. The Path Integral Formulation 

For the numerical simulation we use the path integral formalism, as described 
in Ref. [0. We begin with the Trotter break-up 

e-^^=(e~— j (18) 

and label each of the Nt factors by an imaginary time index t running from 
1 to A^i. One can then split the hopping term from the interaction term 
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making an error of order 



kUJ3^ 



As a next step we introduce a scalar field 



{(Jx,t} for the time slice t in order to eliminate the four fermion term 



e 



^ J2 (4s + 

* {x,y) 



(20) 



Furthermore, we use the two relations 



exp 



l + a\ay,x^y and exp 



l+aK[e"-l] , (21) 



and the same relations for 6^ and b to transform the sums in the exponential 
ofEq. (M) 



exp 



J2 {(^l^i 



:22] 



and 



exp 



I Up 



Nt 

J]^ 1 + al.ax exp 



ax t + const. 



iup_ 



ax t + const. 



(23) 



For every time slice we put the appropriate factors together and use the 
anticommutation relations to transform the expression into a sum of normal 
ordered operators, neglecting terms of higher order in 1/Nt. The trace of 
the product of normal ordered expressions in the fermion operators can be 
transformed into Grassmann integrals in the standard fashion [0], using 



Tr 



■Nt -1 „ Nt 

n : /. {a\ a): = / [] [(dr/*) ft {Vt, Vt) (dr^t) e*^*^''-'''-)] (24) 
.1=1 J t=i 



with the boundary condition rjo = —rjj^^ and similarly for fi{h\ b). Finally, we 
obtain the partition function Z, up to an unimportant numerical coefficient, 
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as 



x,t 



x,t 



x,t 



x,t 



exp 



N't 



iVx,tVy,t + Vy,tVx,. 



{x,y),t 

x,t Vx,t~l J 



x.t 



+ Vx,tVx,t I exp 

x,t 



- 1 



(25) 



when the 6,6 Grassmann variables are associated with 6^,6. This defines a 
fermionic matrix M, specified for each set of cr-fields by its elements between 
the vectors f]y^t' and rjx^t 

J2 iv^,tVy,t + Vy,tVx,t) - J2 V^,tVx,t-l 

(x.v).t x,t 



7]Mr] 



+ Y Vx,tVx,t(iW 



x,t 



(26) 



The integrals over the Grassmann variables may be performed, and we 



get 



YlldcTxA e--^^^.^<^ (detM)' 



(27) 



x,t 



Furthermore detM = (detM)^ so that (detM)^ = det (^M^^mV which is 
useful for the construction of the algorithm. It is possible to include the terms 
of order in the action. They are given in Appendix B and included in 

the simulation. 

In order to compute expectation values of observables, it is convenient to 
first consider that of the operator 



O = exp 



exp [ja] 



{2i 



when j, j are Grassmann sources for the operators a, taken at equal time, 
say t = 0. Inserting O as an additional factor into the l.h.s. of Eq. (p^) 
corresponds to insert the factor O into its r.h.s., where 



O 



driNt+idr}Nt+i exp [jr/jv,+i + Wt+ij] 
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X exp [riNt+i iVNt+i - VnJ] 
X exp[f]i{r]N,+i-VNt)] ■ 



(29) 



This is obtained by associating to [a, a^j an additional Grassmann pair 
{r]Nt+i,VNt+i)y taking care of the boundary conditions which now says rjo = 
—rjNt+i instead of r/o = ~VNt (this explains the last exponential above). 
Explicit integration yields 



O = exp [jT]Nt + Vd + Jj] 



(30) 



In fact O could have been inserted in the trace between any two successive 
time slices t and t + 1 (with the same result for the trace due to cyclic 
invariance). The resulting O would then read 



O = exp [jr]t - rjt+ii + jj] 



(31) 



This technique generalizes to the computation of any (time dependent) 
correlator C,. After C, is rewritten in normal order (any creation operator 
brought to the left of any annihilation operator by repeated use of their 
commutation relations), one makes the replacement in C, d*^^} l*^}^ 



[al[t\,ay [u]] 



d 



d 



'-'Jx,t+i '-'jy,u 



(32) 



the similar one for the 6's, and ends up with the final result: 

^^ id 



(C) = C 



9 j j ' 1 dj 



(33) 



i=j=o 



where 



n 



da^.texp -- 



(T, 



x,t 



X 



exp 



X (detM)^ X 



x,t 



(34) 



For the observables of interest, listed in Eqs. (^, |^ - |^), the following 
expressions then follow. 



in) 



2-2(M-^^^,,i 



P{x ~y) = bxy 



1-2(M-^^^,,i 



(35) 



(36) 
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^^^^ = 2<5.,(m-^^^,,,)-2(m-^^^^,,M-^ , (37) 



C{x -y)= 25.y (M-^^^,,i) + 4-8 (M,;^^^,,,; 

+ 4(m-^^^,,iM-^^;^,i)-2(m-^^^^,iM-^^^,,i) . (38) 

We will also consider the average value of the field a'^ defined as 




Indeed, from the equation of motion for the a-fields, which can easily be 
derived from Eq. (^) or (^), we have 

hm {a',,) = VU{n) . (40) 

This relation, together with (n) = 1 at /i = in the large Nt limit (at finite 
Nt the particle-hole symmetry is broken 0]), will be used as a check that we 
actually considered large enough Nt values. 



4. The numerical simulation 

The expectation values of ( |35| - p9D are obtained from cr-field configurations 



corresponding to the partition function (pTf). These configurations are gen- 
erated with an Hybrid Monte Carlo algorithm (HMC-algorithm) which is a 
widely used tool to simulate systems involving the fermion matrix (e.g. 0, 
0). Its main virtues are that it is an exact algorithm and that, at least 
in principle, the computer time grows slowly with system size. 

Introducing one set of pseudofermion fields 0^*^ for each species r] and 9 
and momenta p conjugate to the a fields, the Hamiltonian of the system can 
be written as 



n 



^Ep? + ^E-l + EE</'?K^);>i' ^ (41) 



where j and k denote sites on the x Nt lattice. Starting from a configu- 
ration a, a new configuration is obtained in three steps. 

i) First the pseudofermion fields are generated from gaussian distributed 
vectors r*-*^ by 0^*^ = M'''((T)r*^*^ and gaussian momenta p are chosen. 
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ii) Then o and p are updated by molecular dynamics with 7i as Hamilto- 
nian, the set of Hamilton's equations being solved with a discrete time 
step dr. So 7i is a constant of motion in the limit dr ^ 0. The inte- 
gration of the molecular dynamics equation is performed by a leapfrog 
method [Bl 



C'"il(n+l)dT = C^jlndr + drpj 



(dr)2 m 



2 do^ 



ndr 



|(n+l)dT 



Vj I uAt 



dr / m 



2 \daj 



dU 



ndr i 



(n+l)dTy 



(42) 
(43) 



iii) After tt-md steps, the resulting a configuration is accepted or refused by 
a Metropohs test. 

This algorithm introduces two parameters, the number of steps timd and 
the step size dr. The correlation between configurations decreases when 
the trajectory length r = nnodr increases, while increasing iimd increases 
the computation time and increasing dr decreases the acceptance. We have 
taken r ~ 1.0, (more precisely between 0.8 and 1.6) with ^md such that a 
large enough acceptance is reached for all the volumes used (the acceptance 
decreases with increasing volume |Q,[§]). This is obtained with tzmd ~ 4/5 [/. 
We also tested trajectories with Poisson-distributed lengths, i.e. dis- 
tributed around a certain mean value and dr held fixed |TD|, O]. This does 
not decorrelate measurements of observables compared to r fixed. The only 
effect is to increase the acceptance a few percent if tt-md is small. 

For an observable O whose value at Monte Carlo time is Oi the autocor- 
relation function is 



We estimated the autocorrelation time by 



tc-l 



o 



Co{t + l) 
Co{t) 



(44) 



(45) 



disregarding the noisy region t > where Co{t) is less than 0.05. This gives 
the exact autocorrelation time if Co is purely exponential. The integrated 
autocorrelation time 



tc-l 



tint — X! ^oit) 



(46) 



coincides with within 5% under the same assumption. 
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becomes less than 10 times the value obtained for the 



The physical and algorithmic parameters are given in Tables (00) along 
with the number of trajectories simulated, the autocorrelation time C,a' and 
the acceptance. Other autocorrelation times connected to fermionic quanti- 
ties are always very small as compared to ^o-'- With the parameters chosen, 
the smallest values of the acceptance at [/ = 4 are 0.64 for the lattice 8^ x 64 
at /? = 6 and 0.71 for the lattice 8^ x 32 at /3 = 4. 

The updates of aj and pj (Eqs. ^) require twice the computation 

of a vector u = for a given (f). This computation, which costs 

most of the computer time, is performed by a conjugate gradient subroutine 
with a diagonal preconditioning [|12], which speeds up the convergence con- 
siderably. The conjugate gradient iteration is stopped after Ccg steps, when 

initial guess Uq. 

The number of conjugate- gradient iterations can be significantly reduced 
if we use for uq the value obtained at the previous molecular dynamics step. 
However this method leads to dramatic decreases of the acceptance when 
increasing f] and V. So we always start with mq = 0, in such a way that no 
peculiar direction of phase space is favored or suppressed by the numerical 
uncertainties. 

The number of conjugate gradient iterations needed to reach the pre- 
scribed accuracy ranges between ~ 30% and ~ 65% of the vector length. 
However in some cases, this accuracy is not obtained in x Nt iterations, 
in particular for small volumes. Such situations can be avoided by reducing 
the trajectory length r, and this explains the values r = 0.8 and 1.2 in Table 
(|l|). However, this comes with the price of increasing the cr autocorrelation 
times. 

The matrix M used in this computation is given in Eq. (|9lD of Appendix 
B. It is valid to order l/N^. Some of the effects of using this second order 
matrix instead of the first order matrix of Eq. ( p6D are shown in Table (^. At 
/i = 0, and for A^^^ large enough, one expects (n) to approach 1 and, according 
to Eqs. (^,1^), (^CTx^t) to approach y/U. This is clearly better achieved with 

the improved action (|9T|) . While Pq is not affected within errors by the change 
of the action, x unexpectedly happens to depend more on Nt in the second 
order case. 

Because increasing Nt is costly in computer time, we currently used Nt = 
2(311. The exceptions, apart for analysis of Nt dependencies, concern the case 
P = 0.5, U = 8 {Nt = 32) and the low temperature points (3 = 6, U = 4 
{Nt = 64). 

The results obtained for (a'), filling of the band (n) /2, pairing correlation 
Pq, spin susceptibility x single occupation probability Si for various 
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lattices and values of the inverse temperature (3 and chemical potential 
are given in Tables The errors reported therein are calculated by 

incorporating the autocorrelation time tint of Eq- (0) by using a variance 
which is the naive variance multiplied by (1 + 2tint). 

We cannot compare the efficiency of our algorithm with Quantum Monte 



Carlo [13 1, because for that algorithm we do not have the information needed. 



5. The S-Wave Pair-Field Correlation Func- 
tion 

One of the observables of interest in the context of a superconducting phase 
transition is the pairing correlation function P and its Fourier transform P. 
They are defined by ([7p and (|T^ respectively. In the present work we have 
studied the q = (0, 0) equal-time s-wave pair-field correlation function Pq, Eq. 
([T3|) , which has previously been examined e.g. in @, ||. For this quantity we 



obtain from Eq.(p6D 

Po = 1 - 2 (M-l,^^ + 2E (K:^.;x,i)') • (47) 



Let us recall |^ what can be expected on general grounds for the 
behavior of Pq as a function of (3, U, fi and A^^, the linear size of the lattice. 
A first remark is that under a global particle-hole transformation 

a. ^ i-Tal (48) 

b. - {-)% (49) 

which changes H{^) of Eq. (|l|) into H{—fi), up to an irrelevant additive 
constant, the operator defining P in Eq. (^ remains unchanged. Hence 

Poifi) = Po{-fi) , (50) 

whereas the filling density verifies 

(n) (^) = 1 - (n) (-/i) . (51) 

In the vicinity of = 0, and above any transition temperature which 
could lead to non-analyticities, we thus expect Pq to be a function of /i^, 
while (n) (fi) is linear in /x. 

At /i = 0, because of the 5'f/(2)' symmetry (Eq. |^), no transition may 
occur at T 7^ 0, and as T ^ 0, the correlation length ^ and Pq are expected 
to diverge according to 

e oc exp (A/T) (52) 
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Po-be ■ (53) 

At yU 7^ 0, the SU{2y symmetry is broken down to a U{1) symmetry 
associated with the particle number, which can lead to a finite T transition 
of the Kosterlitz-Thouless type. Then the expectation is, above the critical 
temperature Tkt, 



^KT oc exp 

and 



A/^T-TKTifi,U) 



(54) 



^^kt""" , (55) 



with r/KT = .25. Of course A and b are functions of fi and U. Also T 
dependent prefactors may occur in these expressions. But it is reasonable 
to think that the fi,U dependences of ^kt and Pq are dominated by that of 
T — Tkt- 

In a finite system of size A'^^, standard scaling arguments predict a be- 
havior 

Po_^^ ~ Ar2/(iV,/0 at ^ = , (56) 

and 

Po,7v. ^ Ar2-''KT/(Ar,/^KT) at ^ ^ (57) 

for the pair correlation, in a domain of parameters (3, fi, U where the ^'s are 
large, and Ng/i is kept fixed. 

On the basis of numerical QMC simulations of the model, Ref. con- 
firmed this overall picture and concluded that a maximum transition tem- 
perature of order 0.2 is reached for f/ ~ 8 and iji) ^ .85, while at f/ = 4 the 
critical temperature Tkt is of order 0.1 and 0.05 respectively at {n) = 0.87 
and 0.5. 

We will now discuss the results of our numerical HMC simulations. We 
start with our results at ^7 = 4, where we have most of the data. 

We first focus on the /^-dependence of Pq, more convenient than its filling 
dependence because (n) (fi) suffers from finite size effects which are not of 
the same nature as those we are interested in, being important even in the 
non interacting case for the Ng values commonly considered. At fixed T and 
U, the above discussion indicates that Pq is fiat close to /i = 0, and should 
increase with \fi\ as a consequence of T — TKT(Ai, U) decreasing, a maximum 
of Tkt leading to a maximum of Pq in the same fi region. This qualitative 
behavior should survive in a finite, large enough box provided the function 
f{x) of Eq. (|F7|), which is proportional to x"^ ''^ at large x, stays monotonic. 

In order to probe such a behavior, we select out from our data of Table (^) 
the values of Pq^Ns obtained at the largest /? (4 and 6) and Ng = 6 available 
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for a sufficient range of fi values. The result is shown in Fig. (|I|). These 
data are compatible with Po,Ns being constant with respect to fi, and not 
compatible with a sharp maximum in this n range (corresponding to a filling 
between .65 and 1.0 for both /3 values), as suggested in Ref. 0. 

We may still examine the dependence of Pq^Ns on at fixed /x, to see if 
these data give evidence for a KT transition at /i 7^ 0. We therefore fix /i at 
-0.15, a value at which we have data in the range 1 < /5 < 6 and for various 
Ng values in the 4-8 range. 

Let us first assume that at each f3 the largest Ng value available is 
large enough for the system to be close to its thermodynamical limit, an 
assumption which is verified at /3 = 4 as seen from our data at A^^ = 4, 6, 8. 
In Fig. (^) we show Po,Ns as a function of T for A^^ = 4 at /3 < 3 and 
Ng = 8 for (3 > 4, together with a fit to Eqs. (0, ^5]) (continuous line). 
The fit is remarkably good in the entire temperature range. It favors a KT 
transition at Tkt = -03(2), at about 4 standard deviations away from the 
value .1 preferred in Ref. Also shown is a fit to Eqs. (|52| , ^3]) (dotted 
line), which is compatible as well with the data, provided the T = 1 point is 
ignored. Imposing Tkt = 0.1 gives poor fits (see Table @). 

These fits provide us with a guess for the parameterization of the corre- 
lation length, up to a normalization factor, which allows us to test scaling 
behavior given by ( |57| ) or (^Bj) . With the parameters of lines 1 and 4 of Table 
d^), we obtain the results shown in Figs. (H) and (|) respectively, together 
with the asymptotic limits c(A^s/Ckt)~^'^^ and c'{Ns/^)~^ expected in the 
infinite volume limit for the scaling functions /kt and / (c and c' are deter- 
mined by the = 3, A^s = 4 point). Since both figures may be considered 
as suggestive of a scaling behavior, the evidence for a Kosterlitz-Thouless 
transition at fi = —0.15 remains poor. 

From this analysis, our conclusion is: If there is a Kosterlitz-Thouless 
transition for U = 4, /i = —0.15, it occurs at such a small value of T in 
units of the hopping parameter k (a few percent) that its effects are hardly 
distinguishable from those of the T = singularity expected at = 0. (Eq. 
^. 

This conclusion is different from that of f^, and one must try to under- 
stand why. At the quantitative level, for U = A, the values we find for Po,Ns 
at /i = —0.15 above /5 = 3 are smaller than in at {n) = 0.87, especially at 
(3 = 6. That (n) depends on A^^ at fixed fi as already mentioned (see Table 
(0) ) is insufficient to explain the difference since we do not observe significant 
fi dependence in -Po,Ar^. 

Of course larger Pq values at a given (3 favor larger Tkt when analyzed in 
terms of the (|57D Ansatz, and this may explain the value Tkt ~ -1 proposed 
in 0. However, we would like to point out that the data corresponding to 
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the highest (3 values 8 and 10 of the latter reference were not used in the 
finite size analysis provided there. A look at these data strongly suggests 
that including them would have considerably lowered the guessed value of 

We did not take enough data at f/ = 2 and [/ = 8 to allow for an analysis 
similar to that achieved at [/ = 4. A look at Table (|]) however confirms that 
the pair correlation increases with /? at fixed [/i, ?7], while its values at /i = 
and -0.15 are consistent with each other at fixed Some information 

comes from comparing Pq at f/ = 2, 4 and 8 for {3 = 2 and /i = 0. fixed. Its 
values are respectively 1.127(4), 1.63(3) and 2.2(2). It would be incorrect to 
directly infer from this that also the correlation length increases. In fact, as 
U increases at fixed temperature, more pairs are formed and this affects the 
normalization of the correlation function. As a measure of this normalization 
we take the x = y (zero distance) value of the correlation of Eq. (|^), 

^(0) = {cilb%a, + b.a^albl) , (58) 

which is nothing but the average number densities of doubly occupied and 
empty sites, related to the single occupation probability Si, given in Tables 
(^,|), by -P(O) = 1 — ^i. Hence a better representative of the correlation 
length variations is Po/P{0), and for this quantity at /? = 2, /i = 0, ?7 = 2, 4, 
and 8 respectively we find 1.80(2), 2.18(4) and 2.4(2). 

This leaves uncertain the guess that the correlation length is larger at 
U = 8 than it is at f/ = 4. 

The above remark about the normalization of Pq, which might be impor- 
tant for a discussion of the U dependence of the KT transition temperature, 
does not alter our previous analysis at ?7 = 4 fixed because, as seen from 
Table (^), -P(O) is essentially insensitive to /3, /x and A^^ in the explored 
ranges. This is true even at temperatures as high as 1 at f/ = 2 and 4. The 
largest variation observed for P{0) is found between T = 0.5 and T = 2, for 
U = 8, where it however does not exceed 10%. So, in the U ( and /i) range 
considered, the pairs are very robust against temperature rises. 

The independence of Si on the filling for fixed U in the region we investi- 
gated means that the increase in the filling goes into the formation of pairs. 
This is, however, also true in the free case in this region. This may be related 
to the fact that 5*1, as well as Pq, is an even function of /x. 
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6. The susceptibility: Comparison of numeri- 
cal results with hopping parameter expansion 

The free energy of the repulsive Hubbard model has been expanded in powers 
of the hopping parameter k by various authors [H, IB, ^. We start from 
the result of 0], where the expansion is pushed to order 5 in /c^. Via a 
particle-hole transformation on one of the two electron operators, the free 
energy expansion is transformed into that relevant for the Hamiltonian in 
the attractive case Eq. (|l]). For arbitrary values of (3~^, /i and h respectively 
of the temperature, chemical potential and external magnetic field, we end 
up with the following truncated expansion ( F is the negative of (3 times the 
free energy density) 



F{f3,U,fi,h) 



1 

V 



\nZ + 0{k 



12\ 



,2k 



log^; + — 
hum' 



(59) 



where 



X 
V 

z 



e^^, w = exp 



1 + xw {y + 1/y) + x'^ 



/3U' 



(60) 



and {R} represents the set of integers {i 



5 ''^J.■> 



The symbol Ck,{_r} represents the 3836 non zero coefficients provided in 
labeled by the order n in k"^ and the associated set {R}- By differen- 
tiation with respect to /i and one obtains similar series (here taken at 
zero magnetic field) for the filling density, the static uniform charge density 
correlation and the static uniform spin susceptibility. We focus on the latter, 
and study the truncated series 



X(g = 0) = [y 



'E -^nk 

K=0 



2k 



(61) 



The coefficients follow from the expansion (^) of F. Although it is 
not strictly speaking a high temperature series, but a series in (/9fc^) with 
/^-dependent coefficients, extracting information from such short expansions 
at low temperature is notoriously difficult. 
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An idea of the problem at hand is given by the following examples. We 
take [/ = 4, /i = and compute the coefficients for (3 = 1 and 2. Their 
numerical ratios to Aq are shown in Table (|^). 

One observes a quite discouraging pattern of orders of magnitude and 
signs (remember we want to evaluate x at = 1). We tried to understand 
what happens by starting from the free case [/ = 0, at /x = 0. The exact 
susceptibility is given by (see appendix A) 

^ I 1 + cosh(2A;/35p) ^ ' 

with 

5p = -(cospi + C0SP2) • (63) 
Its expansion in [jSk)'^ to S*'^ order is found to be 

^ = 1 _ ^pkf + -iPk)' - -{f3kf + ^{(3kf - ^{PkY' (64) 
p ; ^ 2^-^ ^ 36^^ ^ 288 ^'^ ^ 800 ' ^ ' 

The signs now alternate, but again the size of the coefficients in k"^ in- 
creases fast, the more so the temperature is lowered. Here the origin of this 
blow up is clear from expression (^2]). Rewriting it as 

_ 2 
X 



P I ^^^(^)l + cosh(2W 

where p{5) represents the density of states, and recalling that p has only 
a logarithmic (van Hove) singularity at 5 = 0, all singularities of x in ^ 
at finite distance are end point (5 = 2) singularities. They are located at 
cosh(4/c/9) = —1, and the closest one, which fixes the radius of convergence, 
is at {(3kY = — 7r^/16. Hence the alternate signs in the above series and also 
the fact that 7r^/16 times the ratio 1 745/^44! is very close to 1 (namely .990) . 
Finally, we note from the above integral that the most singular part of xl P 
is proportional to tanh(2/5/c)/(/3A;). This leads us to perform the change of 
variable pk ^ u defined by 

u = tanh {2pk) /(pk) - 2 , (66) 

and to re-expand in powers of u. We find 

2x _ ^ 3m 9m2 246^3 4096m^ 30034^^ 

y~ ^ Y ~ 640 ^ 51469 " 2238401 ^ 40581943 ^ " ' ^ '* 
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that is 

^ - 1. + .375M-1.41 10"V+4.78 10~V-1.83 10"V + 7.40 10"V (6 



The interval l3k > is mapped onto — 2 < -u < 0, and low temperatures 
mean u ~ —2. Clearly this series in u looks much more tractable than the 
original one. Indeed, actual comparison of Pade approximants of (|68D with 
the exact result ( ]62| ) shows that, up to /3 ~ 10, the latter can be approached 
with an accuracy better than 5%, an achievement not expected from the 
series (|51D! 

In the absence of any information on the analytic structure of x in ^ 
as the interaction is turned on, we decided to perform the same change of 
variable ([66| ) in all cases and constructed the Pade approximants to the series 
in u. We observe that the [2,3] and [3,2] approximants systematically lead 
to very similar results. Our results for /i = are collected in Fig. (|]), 
and compared with those of the numerical simulation, as obtained with the 
largest lattices which were simulated for each set of (3, U values. We now 
describe and comment this figure. 

The upper (dotted) curve is the exact result in the free case U = 0. The 
four continuous lines are the [2,3] Pade results for respectively U = 0.1, 2, 4, 8, 
from top to bottom . The lower (dashed) line is the "atomic" limit [k = 0) 
drawn for U = 8. The statistical errors are smaller than the size of the 
symbols representing our numerical results for the largest available values of 
Ng, Nt. Tables (|^ |) show that sizable finite size effects are present. In 
particular, unlike Pq, the susceptibility appears to be rather sensitive to the 
number Nf of time slices. When comparisons can be made, our results are 
in general agreement with published data |1^, |^, |T^. 



The U = Q.l curve is used as an overall check that the whole series 



constructed from ||T^ is correctly implemented and that the change of variable 
and Pade reconstruction work close to the free case. For U = 2, excellent 
agreement between series and data is obtained at /3 = 1 and 2. The figure 
suggests that at /3 = 3, the series still gives a reasonable answer, although 
lower than the true one by probably 10 to 15 %. A similar departure shows 
up above j3 = 2 at U = 4, and, say, P = 1.5 at U = 8. Hence the overall 
pattern is that our treatment of the series gives good results at P values 
below a maximum which decreases as U increases. Above this maximum, 
the Pade approximant drops down rapidly the more so when U is large. It 
is interesting to discuss this latter feature. Consider again the case U = 8. 
The dashed line (fc = 0) is the curve (see appendix A) 

I = ,w = exp [-(3U/2] , (69) 
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drawn for U = 8, thus falling down like exp [—4/5]. One observes that the 
Pade result at f/ = 8 drops down roughly parallel to it above (3 = 2. 

Inspection of the series for x reveals that all (known) terms of its ex- 
pansion in k"^ contain at least one power of as a factor. So any truncation 
unavoidably leads to an exponential fall off in j3 at fixed [/, the more so when 
U is large. This is what we observe with the series at hand, while the numer- 
ical simulation indicates that, although x actually decreases as (3 increases, 
it does so at a much slower rate than naively expected. In other words, there 
are collective effects which rather tend to maintain a sizable susceptibility, 
even though low temperature or/and large U favor the formation of local 
pairs, which contribute zero to the total spin. We note that, as already men- 
tioned in the previous section, the probability for single occupation, 5*1 (as 
given in Tables (^^) and thus the average local pair number {n/2) — 6*1/2, 
depends weakly upon temperature. Hence the decrease of x with T is not a 
consequence (only) of pair formation at fixed U and (n). Its sensitivity to 
U also is strong compared to that of 5*1. As recently discussed in [|18[, these 



questions are of interest in connection with the kind of regime (BCS-like or 
pair condensation) eventually leading to a superconducting transition. See 



also [H, 17 



Our results show that, unlike popular approximations such as RPA or 
T-matrix treatment, the series expansion in k"^ can be continued successfully 
down to fairly small T values even at large coupling. In order to further 
illustrate this point, we present our predictions for x as a function of T for 
U = 4, 6, 8 and 12, at half filling ((n) = 1), in Fig. (6). These predictions 
appear as solid lines when our calculation is accurate (the minimum T's of 
the solid lines are conservative guesses). We point out that by themselves 
these results clearly expose the main features of physical interest, namely the 
existence of a crossover from a smooth regime at high T to a rapid fall off 
below a U dependent T value. Furthermore, we observe that these x values 
are quantitatively very close to that of |]18| obtained numerically at the same 
U values, but for {n) = 0.8. This enforces our general statement that there 
is no evidence for any substantial dependence on filling. 



7. Summary and conclusion 

Motivated by the interest for strongly interacting systems of electrons in 
the context of high Tc materials, we have performed an investigation of the 
attractive Hubbard model. 

Based on the hybrid Monte-Carlo method to simulate the path integral 
representation of the partition function, a numerical study provided us with 
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estimates of the static zero-momentum pair correlation Pq and spin suscepti- 
bility X- The exploration covered the ranges [0.5, 6.0] for the inverse temper- 
ature /3 and [2.0, 8.0] for the interaction strength U in units of the hopping 
parameter k. At 11=4: we investigate in detail the dependence on the chemi- 
cal potential in order to study the sensitivity of these physical quantities to 
band filling, at and away from half filling. We also performed an analytical 
investigation of the susceptibility, based on its expansion in k'^ to order 5 
derived from the results of 0] and on a specific Pade resummation technique 
inspired by the free case. We show that this method yields very good results 
for /3 < 1 at any U, and /5 < 2 for U < A. 

In the numerical approach, the cost in computer time leads to the con- 
sideration of rather small lattices only, and also limits the number of time 
slices involved. The limitation of the analytical approach is obviously due 
the truncation of the fc^ expansion, which furthermore is not available for the 
pair correlation. But it yields estimates of x directly in the thermodynamical 
limit, so that the two approaches are complementary and their comparison 
brings valuable information. Whenever it was possible, we also compared 
our numerical data with those of previous investigations |l8l performed 
using Quantum Monte-Carlo. The overall agreement is good, although we 
found significant discrepancies for Pq our lowest T values. 

Using pair correlations, the main question addressed was that of the exis- 
tence and location Tkt of a Kosterlitz-Thouless superconducting transition, 
the only possible transition at finite temperature away from half filling. Our 
results show very weak dependence on /i, hardly distinguishing half filling 
from other values of band filling in the range between .4 and .5. Fits or finite 
size analysis relevant to a KT transition are consistent with a transition at 
Tkt/^ of the order of a few percent, with a weak dependence on fi, if any. 
This conclusion is at variance with that drawn in but we point out that 
it is perhaps consistent with what might follow from taking into account the 
lowest temperature data points of this latter reference. In any case, from our 
analysis we conclude that it is always necessary to compare quantitatively 
the behavior at yU = 0, where there should not be a transition at finite T to 
fi 0, before one can ascertain the existence of a KT transition at finite T, 
and estimate its critical temperature. 

Our study of the magnetic susceptibility shows that it depends very 
weakly on the band filling, but decreases rapidly with increasing U and/or 
decreasing temperature, in agreement with previous investigations ^7\, ||, |18[. 
Increasing U favors S-wave pair formation, reducing the probability 5*1 of sites 
occupied by a single electron responding to a magnetic field. This explains 
the simultaneous decrease of x and S'l at fixed temperature. But at fixed 
U, we observe that 5*1 stays approximately constant in the low temperature 
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range whereas x faUs steeply. This casts some doubt on the existence of a 
tight connection between susceptibihty and pair formation as the expected 
transition is approached as claimed in ||18[ |. 

In this work, we restricted our data analysis to static, zero momentum 
correlations. However, we have stored all the information required to also an- 
alyze spatial correlations and investigate dynamical properties of the model. 



Acknowledgments 

We would like to thank Thierry Jolicoeur for many interesting discussions 
and for his careful reading of our manuscript. We further thank J. Oitmaa 
for providing us with the latest results on the series expansion in suitable 
form. Two of us (B. P. and J. S.) are grateful to the Deutsche Forschungsge- 
meinschaft for support. Calculations were performed at HLRZ (Jiilich) and 
CEA (Grenoble). 



Appendix A 

We first discuss the atomic case k = 0, at h = 

H = -U (a^a - (^b^b - ^) - (a^a + b^b) . (70) 

There are only 4 states. It is convenient to introduce the variables 

X = e^'^ , (71) 

w = e-— . (72) 



One easily finds 



Z = l + 2wx + x'^ , (73) 

2x {w + x) 
1 + 2wx + 



. (74) 



f.- , , , (75) 



1 + 2wx + x^ 

^ = ^ l + 2wx + x^ ' ^^'^ 

C(0,0)= 2x(^ + 2x) _ 
^ ' ^ l + 2wx + x2 ^ ^ 
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The single occupation probability is 



2xw X 
' " 1 + 2WX + " ■ ^ 

One remarks that Pq is always smaller than 1. The path integral expres- 
sion of Z factorises into one path integral for every site. It can be analytically 
calculated, and with the form of the action which we use, the result is exact 
independent of the value of Nt, i.e. there are no finite Nt corrections. 

In the free case U — 0, the Hamiltonian, diagonalized in momentum 
space, is 

Hu=o = ^{ep - n) {dldp + - /i(ajap - 6j6p) , (79) 
p 

where = — 2A;(cospi + C0SP2) for momenta p = (pi,P2) and h is the mag- 
netic field introduced for completeness. 

On a finite lattice pi = j^rii, rii — 1,2, • ■ • , Ng. The correspondence with 
the thermodynamical limit Ng — > 00 is given by 

N^^ Jo 47r2 ^ ^ 

Setting 



Sp = exp[-/5(ep - fi)] , (81) 
one derives the partition function 

p 



and the zero field quantities of interest : 



1 ^ 1 + s 



2 



p(l + 2s, 

{l + Spf 
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The single occupation probability 5*1 = {n) — {n)'^ /2 is independent of 
temperature at fixed filling. It varies from at zero filling to 1/2 at half 
filling. 

Using Sp = l/sp_(7r,7r) at /i = 0, the susceptibility at half-filling becomes 



fj.=0 



m ? 1 



s p 



+ cosh(/5ep) 



^7) 



Appendix B 

In order to include terms of order within the action, we apply a Trotter- 
splitting in the following way: 



where 



Hk = {(^IcLy + clIcI'x + blby + h^J^ 



and 



(89) 



(90) 



We then, as before, introduce an auxiliary scalar field. Then we expand the 
three exponentials in (|88|), collect terms including order and write the 

summands of operators in normal order. In the resulting expression we can 
substitute the operators by Grassmann variables. The procedure described 
gives rise to additional terms in the action. The spatial elements of the 
matrix in Eq. ( p6D now read 



y,t 



k(3 
2N, 



J2 i^y,^+i> + 



y,x—u J 



exp 



+ 5^^y I exp 



+ exp 



+ 



^ 7/ !>' 



(91) 



i>,u'=l 
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Table 1: Parameters of the simulations performed at [/ = 4. Also given 
are the number of trajectories generated, the a' field autocorrelation time ^o-/ 
and the acceptance. 
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Table 2: Parameters of the simulations performed aX U = 2 and 8. Also 
given are the number of trajectories generated, the a' field autocorrelation 
time and the acceptance. 
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Table 3: Comparison of data taken with the fermionic matrix to order 1 and 
2 in l/Nt- See Eqs. ( p6D and ( [91| ) respectively. 
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42 X 16 


4 


2 


.0 


2.02(3) 


.494(4) 


1.63(3) 


1.18(2) 


.244(2) 


42 X 24 


4 


3 


.15 


1.81(3) 


.435(3) 


2.18(5) 


0.85(3) 


.234(2) 


42 X 24 


4 


3 


.0 


2.10(3) 


.509(6) 


2.23(5) 


0.87(2) 


.242(2) 


42 X 32 


4 


4 


.6 


1.31(2) 


.309(1) 


2 47(7) 


0.45(2) 


.226(4) 


42 X 32 


4 


4 


.45 


1.40(3) 


.335(2) 


2 6(1) 


0.46(2) 


.234(1) 


42 X 32 


4 


4 


.15 


1.87(6) 


45(11 


2.40(6) 


0.55(2) 


.234(4) 


42 X 32 


4 


4 


.0 


2.11(4) 


.514(4) 


2.37(9) 


0.51(2) 


.230(4) 


6^ X 32 


4 


4 


.45 


1.38(1) 


.329(1) 


3 0(1) 


0.72(1) 


.226(1) 


6^ X 32 


4 


4 


.2 


1.74(2) 


.422(2) 


2.87(5) 


0.697(9) 


.234(1) 


6^ X 32 


4 


4 


.15 


1.83(2) 


.443(2) 


3 1(1) 


0.68(1) 


.232(2) 


6^ X 32 


4 


4 


.1 


1.92(1) 


.468(2) 


3 1(2) 


0.71(1) 


.234(1) 


6^ X 32 


4 


4 


.0 


2.06(2) 


.501(3) 


3.05(9) 


0.75(1) 


.239(1) 


8^ X 32 


4 


4 


.15 


1.86(2) 


.449(2) 


2.85(8) 


0.72(1) 




42 X 64 


4 


6 


.45 


1.31(3) 


.321(2) 


2.38(8) 


0.14(2) 


.246(6) 


42 X 64 


4 


6 


.3 


1.41(4) 


.344(2) 


2.9(1) 


0.22(2) 


.242(6) 


42 X 64 


4 


6 


.15 


1.61(4) 


.396(3) 


3.3(2) 


0.23(2) 


.234(4) 


42 X 64 


4 


6 


.08 


1.86(6) 


.460(7) 


2.58(9) 


0.18(1) 


.222(4) 


42 X 64 


4 


6 


.0 


2.03(4) 


.503(7) 


2.57(9) 


0.17(2) 


.236(6) 


62 X 64 


4 


6 


.45 


1.32(2) 


.322(2) 


4.0(3) 


0.38(3) 


.224(6) 


62 X 64 


4 


6 


.3 


1.46(2) 


.359(1) 


4.2(2) 


0.37(4) 


.240(4) 


62 X 64 


4 


6 


.23 


1.62(2) 


.396(1) 


4.2(2) 


0.42(2) 


.244(4) 


62 X 64 


4 


6 


.15 


1.83(5) 


.452(2) 


3.9(3) 


0.31(2) 


.238(4) 


62 X 64 


4 


6 


.0 


2.01(2) 


.493(2) 


3.8(2) 


0.36(3) 


.236(4) 


82 X 64 


4 


6 


.15 


1.75(3) 


.429(1) 


4.9(4) 


0.49(3) 


.242(2) 



Table 4: Data obtained on various lattices at C/ = 4 : {a'), filling of the 
band (n) /2, pairing correlation Pq, spin susceptibility X ^^^^ single occupa- 
tion probability 5*1 for different values of the inverse temperature and of the 
chemical potential. 
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Lattice 


U 


P 


-/i 


{(^) 


{n) /2 


Po 


X X 10 


Si 


4x8 


2 


-I 

1 


.U 


1.45(1) 


.502(2) 


0.886(1) 


2.114(4) 


.3738(8) 


4 X ib 


o 
Z 


o 
/ 


.U 


i.4z(zj 


.49 / (o) 


1 1 07/ /I \ 

1.12/(4) 


2.46(1 ) 


.3/3(1 ) 


6^ X 32 


2 


4 


.15 


1.268(7) 


.438(1) 


1.46(1) 


2.55(1) 


.3732(8) 


6^ X 32 


2 


4 


.0 


1.453(4) 


.504(1) 


1.445(7) 


2.727(9) 


.3746(4) 


6^ X 24 


2 


6 


.3 


1.153(4) 


.3747(8) 


1.99(2) 


1.30(2) 


.3446(6) 


6^ X 24 


2 


6 


.15 


1.340(8) 


.440(2) 


2.20(6) 


1.69(2) 


.3444(8) 


6^ X 24 


2 


6 


.0 


1.605(6) 


.531(2) 


2.10(3) 


1.86(2) 


.346(1) 


6^ X 32 


8 


.5 


.15 


2.69(7) 


.475(8) 


0.948(3) 


0.667(5) 


.174(2) 


6^ X 32 


8 


.5 


.0 


2.88(6) 


.509(8) 


0.955(4) 


0.661(5) 


.173(2) 


6^ X 32 


8 


2. 


.15 


2.3(2) 


.41(3) 


1.9(1) 


0.104(5) 


.079(2) 


6^ X 32 


8 


2. 


.0 


2.7(1) 


.47(2) 


2.2(2) 


0.116(5) 


.085(2) 



Table 5: Data obtained on various lattices at f/ = 2 and 8. As in Table (p. 



fitting 


; function 


b 


a 


Tkt 


Npts 


x'/d.o.f. 


6exp 


a/ a/T - Tkt 




0.42(2) 


0.92(6) 


0.03(2) 


5 


1.58 






0.46(9) 


0.9(2) 


0.04(3) 


4 


2.96 


bexp [a/T] 


0.768(8) 


0.350(7) 




5 


10.55 






0.90(3) 


0.30(1) 




4 


1.73 


bexp 


a/VT-0.1 




0.545(9) 


0.66(1) 


0.1 


5 


16.98 






0.68(3) 


0.55(2) 


0.1 


4 


6.02 



Table 6: Fits of Pq as a function of T for U = A and n = —0.15. The point 
T = 1 is ignored in every second line. In the last two lines, Tkt=0.1 [5] is 
imposed. 



K 





1 


2 


3 


4 


5 


(3 = 1 


1. 


.580 


-.725 


.372 


.156 


-.497 


P = 2 


1. 


5.65 


2.91 


-17.6 


-.404 


62.5 



Table 7: The ratios A^^/Aq of the expansion of x iii powers of k'^ for U = 4, 
;x = 0, /? = 1 and 2. 
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Lattice 


U 


/5 






{n)/2 


Po 


XX 10 


^1 


62 X 32 


2 


4 


.15 


1.268(7) 


.438(1) 


1.46(1) 


2.55(1) 


.3732(8) 


62 X 16 


2 


4 


.15 


1.377(4) 


.4514(8) 


1.67(1) 


2.114(9) 


.3510(6) 


62 X 32 


2 


4 


.0 


1.453(4) 


.504(1) 


1.445(7) 


2.727(9) 


.3746(4) 


62 X 16 


2 


4 


.0 


1.583(4) 


.523(1) 


1.65(1) 


2.212(9) 


.3532(6) 


42 X 16 


4 


1 


.15 


1.81(3) 


.453(5) 


1.062(5) 


1.379(9) 


.264(2) 


42 X 8 


4 


1 


.15 


1.90(2) 


.463(4) 


1.072(4) 


1.289(8) 


.250(2) 


42 X 16 


4 


1 


.0 


1.96(3) 


.490(5) 


1.062(5) 


1.390(9) 


.266(2) 


42 X 8 


4 


1 


.0 


2.05(2) 


.501(4) 


1.071(4) 


1.301(8) 


.252(2) 


62 X 16 


4 


1 


.0 


2.02(2) 


.501(3) 


1.072(7) 


1.382(8) 


.267(1) 


62 X 8 


4 


1 


.0 


2.03(2) 


.496(3) 


1.075(6) 


1.313(8) 


.255(1) 


42 X 32 


4 


4 


.0 


2.11(4) 


.514(4) 


2.37(9) 


0.51(2) 


.230(4) 


42 X 64 


4 


4 


.0 


1.98(3) 


.495(5) 


2.55(9) 


0.71(2) 


.244(2) 



Table 8: Comparison of data taken at different Nt values. 
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Figure 1: Pq as a, function of // for C/ — A. 
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T=1/p 



Figure 2: Pq function of the temperature for U 
The sohd curve is the fit to Pq = feexp a/{T — TktY^"^ 



= 4 and /x = —0.15. 
the dashed curve is 



the fit to Po 
I- 



6exp [a/T]. The parameters are given in hne 1 and 4 of Table 
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Figure 5: x//3 as a function of (5 for different values of the coupling U at 
// = 0. See description in the text. 
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Figure 6: % as a function of T for different values of the coupling U at = 0. 
Prom top to bottom U = 4, 6, 8 and 12. See description in the text. 
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